********************************************************************************
*DYNAMIC IMPACTS OF PRICING GROUNDWATER
*Bruno, Jessoe, Hanemann in JAERE
****NOTE: Summer event study uses different dataset CALLED BELOW*** 
********************************************************************************

clear all
capture log close
set more off

*SELECT OUTPUT DATE
global outputdate = "20231129"		 

*SET DIRECTORY
cd  "D:\Ellen\Dropbox\Pajaro_AgInnovation" 

use "Data\Parcel_clean_yearrun_20230626.dta", clear

log using "Submission\JAERE\Replication_Code\Log\Event_study_log.log", replace

***Robustness ES without Q3 uses different dataset -- rename year_short to year_run and rerun the first event study, change name
*use "Data\Parcel_clean_yearshort_20230626.dta", clear
*rename year_short year_run

*Balance the panel
bysort parcelnum: gen ntime = [_N]
sum ntime
keep if ntime==`r(max)'

*To get total acreage number for counterfactual analysis = 24,419
preserve
collapse (sum) parcel_size_fixed area_acres, by(year)
sum  parcel_size_fixed area_acres
restore

********************************************************************************
*1.1 Generate dummy variables
********************************************************************************

*Generate event periods 
	*year_run=6 is Q42010-Q32011
gen event_time = 0
sum year_run
forvalues j = 1/`r(max)' {
replace event_time = -`j' if year_run == 6 - `j'
replace event_time = `j' if year_run == 6 + `j'
}
sum year_run

*Generate 5 pre period interaction terms
forvalues j = 1/`r(max)' {
gen minus`j' = 0
replace minus`j' = 1 if inside == 1 & event_time == -`j'
}

*Generate 6 post period interaction terms
forvalues j = 1/`r(max)' {
gen plus`j' = 0
replace plus`j' = 1 if inside == 1 & event_time == `j'
}

*Generate treatment event interaction
**The date of the price split = 2010 
generate event_treatment = 0
replace event_treatment  = 1 if inside == 1 & event_time==0


*Drop missing variables
foreach var of varlist plus* minus* {
summ `var'
   if `r(sum)' == 0  {
            drop `var'
           }
       else { 
               }
         }

gen del_parcel = delivered/100
gen total_water = extraction + del_parcel

********************************************************************************
*Figure 4 Panel A: MAIN UNCONDTIIONAL EVENT STUDY DROPPING 2010
********************************************************************************
drop if year_run==5

preserve
xtreg extraction plus* i.year_run#i.county event_treatment i.inside minus5 minus4 minus3 if year_run<12, fe cluster(parcelnum)

*TABLE 6*
esttab using "Tables\Eventstudy_annual_$outputdate.tex", label replace ///
	se star(* .10 ** .05 *** .01) ///
	addnote(Table reports results from event study regression. Standard errors are clustered at the parcel level.) ///
	cells(b(star fmt(2)) se(par fmt(2)))
eststo clear

*WALD TESTS
putexcel set "Tables\pairwise_wald_annual_$outputdate.xlsx", replace
putexcel B1=2011
putexcel C1=2012
putexcel D1=2013
putexcel E1=2014
putexcel F1=2015
putexcel A2=2011
putexcel A3=2012
putexcel A4=2013
putexcel A5=2014
putexcel A6=2015
*2011 
lincom _b[event_treat]-_b[plus1]
putexcel C2=(r(t)) 
lincom _b[event_treat]-_b[plus2]
putexcel D2=(r(t))
lincom _b[event_treat]-_b[plus3]
putexcel E2=(r(t)) 
lincom _b[event_treat]-_b[plus4]
putexcel F2=(r(t)) 
*2012
lincom _b[plus1]-_b[plus2]
putexcel D3=(r(t)) 
lincom _b[plus1]-_b[plus3]
putexcel E3=(r(t))
lincom _b[plus1]-_b[plus4]
putexcel F3=(r(t)) 
*2013
lincom _b[plus2]-_b[plus3]
putexcel E4=(r(t))
lincom _b[plus2]-_b[plus4]
putexcel F4=(r(t)) 
*2014
lincom _b[plus3]-_b[plus4]
putexcel F5=(r(t)) 

mat coeff = [ _b[minus5]\ _b[minus4] \ _b[minus3] \ 0 \ _b[event_treatment] ///
			\ _b[plus1] \ _b[plus2] \ _b[plus3] \ _b[plus4] \_b[plus5]  \_b[plus6] \_b[plus7] \_b[plus8] \_b[plus9]]
		
mat se = [ _se[minus5] \ _se[minus4] \ _se[minus3] \ 0  \ _se[event_treatment] ///
\ _se[plus1] \ _se[plus2] \ _se[plus3] \ _se[plus4] \_se[plus5] \_se[plus6] \_se[plus7] \_se[plus8] \_se[plus9]]

mat upper=coeff+((1.645)*se)
mat lower=coeff-((1.645)*se)
mat x=(-5\-4\-3\-2\0\1\2\3\4\5\6\7\8\9)
mat data=[coeff,upper,lower,x]
svmat data
rename data1 coeff
rename data2 upper
rename data3 lower
rename data4 x


matrix colnames data=beta highCI loCI year
matrix list data
putexcel set "Tables\water_eventstudy_cy_90_$outputdate.xlsx", sheet(annual) replace
putexcel A1 = matrix(data)
clear
svmat data, names(col)

replace year = 2006 if year ==-5
replace year = 2007 if year ==-4
replace year = 2008 if year ==-3
replace year = 2009 if year ==-2
replace year = 2010 if year ==-1
replace year = 2011 if year ==0
replace year = 2012 if year ==1
replace year = 2013 if year ==2
replace year = 2014 if year ==3
replace year = 2015 if year ==4
replace year = 2016 if year ==5
replace year = 2017 if year ==6
replace year = 2018 if year ==7
replace year = 2019 if year ==8
replace year = 2020 if year ==9

twoway (rcap loCI highCI year if year <2016) (scatter beta year if year <2016, msymbol(triangle) msize(small)), ///
ytitle(Extraction (AF))  ///
 scheme(s1color) legend(label (1 90% CI from t-test) label (2 Difference Inside and Outside)) ///
	xtitle(Year) title(A. Annual)  yline(0, lcolor(black)) xlab(2006(2)2014) ylab(-60(20)40)

graph save "Figures\ES_ann_cy_90_$outputdate.gph", replace
graph export "Figures\ES_ann_cy_90_$outputdate.png", replace
restore


********************************************************************************
*WEIGHTED REGRESSION WITH SCALED DEPENDENT VARIABLE
********************************************************************************
{
drop if year_run==5
*UNCONDITIONAL EVENT STUDY
*MAKE SURE 2010 IS DROPPED BEFORE RUNNING THIS SECTION
preserve
xtreg extract_acre plus* i.year_run#i.county event_treatment i.inside minus5 minus4 minus3 minus2 [aweight = area_acres] if year_run<11, fe cluster(parcelnum)

mat coeff = [ _b[minus5]\ _b[minus4] \ _b[minus3] \ _b[minus2]  \ _b[event_treatment] ///
			\ _b[plus1] \ _b[plus2] \ _b[plus3] \ _b[plus4] \_b[plus5]  \_b[plus6] \_b[plus7] \_b[plus8] \_b[plus9]]
		
mat se = [ _se[minus5] \ _se[minus4] \ _se[minus3] \ _se[minus2]  \ _se[event_treatment] ///
\ _se[plus1] \ _se[plus2] \ _se[plus3] \ _se[plus4] \_se[plus5] \_se[plus6] \_se[plus7] \_se[plus8] \_se[plus9]]

mat upper=coeff+((1.645)*se)
mat lower=coeff-((1.645)*se)
mat x=(-5\-4\-3\-2\0\1\2\3\4\5\6\7\8\9)
mat data=[coeff,upper,lower,x]
svmat data
rename data1 coeff
rename data2 upper
rename data3 lower
rename data4 x

*Produces Column (1) of Table 7; multiply by total acreage to get column (2)
matrix colnames data=beta highCI loCI year
matrix list data
putexcel set "Tables\policy_water_eventstudy_cy_95_$outputdate.xlsx", sheet(annual) replace
putexcel A1 = matrix(data)
clear
svmat data, names(col)


replace year = 2006 if year ==-5
replace year = 2007 if year ==-4
replace year = 2008 if year ==-3
replace year = 2009 if year ==-2
replace year = 2010 if year ==-1
replace year = 2011 if year ==0
replace year = 2012 if year ==1
replace year = 2013 if year ==2
replace year = 2014 if year ==3
replace year = 2015 if year ==4
replace year = 2016 if year ==5
replace year = 2017 if year ==6
replace year = 2018 if year ==7
replace year = 2019 if year ==8
replace year = 2020 if year ==9


twoway (rcap loCI highCI year if year <2016) (scatter beta year if year <2016, msymbol(triangle) msize(small)), ///
ytitle(Extraction - weighted (AF/acre)) ///
 scheme(s1color) legend(label (1 90% CI from t-test) label (2 Difference Inside and Outside)) ///
	xtitle(Year) yline(0, lcolor(black)) xlab(2006(2)2015) 
	
graph save "Figures\ES_weight_$outputdate.gph", replace
graph export "Figures\ES_weight_$outputdate.png", replace
restore

}




********************************************************************************
*2. SUMMER EVENT STUDY
********************************************************************************
use "Data\Parcel_Clean_20230626.dta", clear
keep if quarter==3

*Balance the panel
bysort parcelnum: gen ntime = [_N]
sum ntime
keep if ntime==`r(max)'


********************************************************************************
*2.1 Generate dummy variables
********************************************************************************

*Generate event periods 
	*year_run=6 is Q42010-Q32011
gen event_time = 0
sum ntime
forvalues j = 1/`r(max)' {
replace event_time = -`j' if year== 2011 - `j'
replace event_time = `j' if year == 2011 + `j'
}
sum ntime

*Generate 5 pre period interaction terms
forvalues j = 1/`r(max)' {
gen minus`j' = 0
replace minus`j' = 1 if inside == 1 & event_time == -`j'
}

*Generate 6 post period interaction terms
forvalues j = 1/`r(max)' {
gen plus`j' = 0
replace plus`j' = 1 if inside == 1 & event_time == `j'
}

*Generate treatment event interaction
**The date of the price split = 2010 
generate event_treatment = 0
replace event_treatment  = 1 if inside == 1 & event_time==0


*Drop missing variables
foreach var of varlist plus* minus* {
summ `var'
   if `r(sum)' == 0  {
            drop `var'
           }
       else { 
               }
         }

drop if year==2010
*Drop Q42009-Q32010 because of anticipation

*UNCONDITIONAL SUMMER EVENT STUDY
********************************************************************************
*Figure 4 Panel B: MAIN UNCONDTIIONAL [SUMMER] EVENT STUDY DROPPING 2010
********************************************************************************

preserve
xtreg extraction plus* i.year#i.county event_treatment i.inside minus6 minus5 minus4 minus3 if year <2017, fe cluster(parcelnum)
*xtreg extraction plus* i.year#i.county event_treatment i.inside minus6 minus5 minus4 minus3, fe cluster(parcelnum)

*TABLE 6*
esttab using "Tables\Eventstudy_summer_$outputdate.tex", label replace ///
	se star(* .10 ** .05 *** .01) ///
	addnote(Table reports results from event study regression. Standard errors are clustered at the parcel level.) ///
	cells(b(star fmt(2)) se(par fmt(2)))
eststo clear


*WALD TESTS
putexcel set "Tables\pairwise_wald_summer_$outputdate.xlsx", replace
putexcel B1=2011
putexcel C1=2012
putexcel D1=2013
putexcel E1=2014
putexcel F1=2015
putexcel A2=2011
putexcel A3=2012
putexcel A4=2013
putexcel A5=2014
putexcel A6=2015
*2011 
lincom _b[event_treat]-_b[plus1]
putexcel B3=(r(t)) 
lincom _b[event_treat]-_b[plus2]
putexcel B4=(r(t))
lincom _b[event_treat]-_b[plus3]
putexcel B5=(r(t)) 
lincom _b[event_treat]-_b[plus4]
putexcel B6=(r(t)) 
*2012
lincom _b[plus1]-_b[plus2]
putexcel C4=(r(t)) 
lincom _b[plus1]-_b[plus3]
putexcel C5=(r(t))
lincom _b[plus1]-_b[plus4]
putexcel C6=(r(t)) 
*2013
lincom _b[plus2]-_b[plus3]
putexcel D5=(r(t))
lincom _b[plus2]-_b[plus4]
putexcel D6=(r(t)) 
*2014
lincom _b[plus3]-_b[plus4]
putexcel E6=(r(t)) 


mat coeff = [ _b[minus6] \_b[minus5] \ _b[minus4] \ _b[minus3] \ 0 \ _b[event_treatment] ///
			\ _b[plus1] \ _b[plus2] \ _b[plus3] \ _b[plus4] \_b[plus5]  \_b[plus6] \_b[plus7] \_b[plus8] \_b[plus9]]
		
mat se = [ _se[minus6]  \_se[minus5] \ _se[minus4] \ _se[minus3] \ 0 \ _se[event_treatment] ///
\ _se[plus1] \ _se[plus2] \ _se[plus3] \ _se[plus4] \_se[plus5] \_se[plus6] \_se[plus7] \_se[plus8] \_se[plus9]]


mat upper=coeff+((1.645)*se)
mat lower=coeff-((1.645)*se)

*mat upper=coeff+((1.96)*se)
*mat lower=coeff-((1.96)*se)
mat x=(-6\-5\-4\-3\-2\0\1\2\3\4\5\6\7\8\9)
mat data=[coeff,upper,lower,x]
svmat data
rename data1 coeff
rename data2 upper
rename data3 lower
rename data4 x

matrix colnames data=beta loCI highCI year
matrix list data
putexcel set "Tables\water_eventstudy_cy_90_$outputdate.xlsx", sheet(summer) modify
putexcel A1 = matrix(data)
clear
svmat data, names(col)

replace year = 2005 if year ==-6
replace year = 2006 if year ==-5
replace year = 2007 if year ==-4
replace year = 2008 if year ==-3
replace year = 2009 if year ==-2
replace year = 2010 if year ==-1
replace year = 2011 if year ==0
replace year = 2012 if year ==1
replace year = 2013 if year ==2
replace year = 2014 if year ==3
replace year = 2015 if year ==4
replace year = 2016 if year ==5
replace year = 2017 if year ==6
replace year = 2018 if year ==7
replace year = 2019 if year ==8
replace year = 2020 if year ==9

twoway (rcap loCI highCI year if year <2016) (scatter beta year if year <2016, msymbol(triangle) msize(small)), ///
ytitle(Extraction (AF)) ///
 scheme(s1color) legend(label (1 90% CI from t-test) label (2 Difference Inside and Outside)) ///
	xtitle(Year) title(B. Summer) yline(0, lcolor(black)) xlab(2005(2)2015) ylab(-30(10)20)
graph save "Figures\ES_su_cy_90_$outputdate.gph", replace
graph export "Figures\ES_su_cy_90_$outputdate.png", replace
restore


*Combine main figs to one 2-panel figure:
graph combine "Figures\ES_ann_cy_90_$outputdate.gph"  "Figures\ES_su_cy_90_$outputdate.gph" , xcommon rows(2) cols(1) 
graph display, ysize(18) xsize(15) scheme(s1color) 
graph save "Figures\ES_combine_$outputdate.gph", replace
graph export "Figures\ES_combine_$outputdate.png", width(4000) replace

log close
